PHYSICS  OF  PLASMAS  13,  023503  (2006) 


Force  interaction  of  high  pressure  glow  discharge  with  fluid  flow 
for  active  separation  control 

Subrata  Roya) 

Computational  Plasma  Dynamics  Laboratory,  Mechanical  Engineering,  Kettering  University, 

Flint,  Michigan  48504 

Datta  V.  Gaitonde 

Computational  Sciences  Branch,  Air  Vehicles  Directorate,  Air  Force  Research  Laboratory, 

Wright  Patterson  AFB,  Ohio  45433 

(Received  12  October  2005;  accepted  19  December  2005;  published  online  24  February  2006) 

Radio  frequency  based  discharges  at  atmospheric  pressures  are  the  focus  of  increased  interest  in 
aerodynamics  because  of  the  wide  range  of  potential  applications  including,  specifically,  actuation 
in  flows  at  moderate  speeds.  Recent  literature  describing  promising  experimental  observations, 
especially  on  separation  control,  has  spurred  efforts  in  the  development  of  parallel  theoretical 
modeling  to  lift  limitations  in  the  current  understanding  of  the  actuation  mechanism.  The  present 
effort  demonstrates  higher  fidelity  first-principle  models  in  a  multidimensional  finite-element 
framework  to  predict  surface  discharge-induced  momentum  exchange.  The  complete  problem  of  a 
dielectric  barrier  discharge  at  high  pressure  with  axially  displaced  electrodes  is  simulated  in  a 
self-consistent  manner.  Model  predictions  for  charge  densities,  the  electric  held,  and  gas  velocity 
distributions  are  shown  to  mimic  trends  reported  in  the  experimental  literature.  Results  show  that  a 
residual  of  electrons  remains  deposited  on  the  dielectric  surface  downstream  of  the  exposed 
powered  electrode  for  the  entire  duration  of  the  cycle  and  causes  a  net  electric  force  in  the  direction 
from  the  electrode  to  the  downstream  surface.  For  the  first  time,  results  document  the  mitigation 
process  of  a  separation  bubble  formed  due  to  how  past  a  hat  plate  inclined  at  12°  angle  of  attack. 
This  effort  sets  the  basis  for  extending  the  formulation  further  to  include  polyphase  power  input  in 
multidimensional  settings,  and  to  apply  the  simulation  method  to  hows  past  common  aerodynamic 
conhgurations.  ©  2006  American  Institute  of  Physics.  [DOI:  10.1063/1.2168404] 


I.  INTRODUCTION 

Various  experimental  efforts 1-6  have  demonstrated  the 
feasibility  of  utilizing  dielectric  barrier  discharge  (DBD)  at 
high  pressure  to  inhibit  phenomena  such  as  stall  on  airfoils, 
thus  improving  their  lift-drag  curves.  A  particularly  attractive 
feature  of  these  plasma-based  devices  is  their  capacity  to 
operate  at  atmospheric  conditions  without  uncontrolled  mac¬ 
roscopic  breakdown.  Experimental  efforts  have  identihed 
many  of  the  key  parameters,  most  prominent  among  which 
are  features  of  the  geometric  conhguration  and  the  form  of 
the  applied  excitation. 

A  schematic  of  a  typical  actuator  is  shown  in  Fig.  1.  Two 
electrodes  are  employed:  the  hrst  exposed  to  the  how  and  the 
second  embedded  in  the  dielectric  and  displaced  in  the 
streamwise  direction  relative  to  the  exposed  electrode.  The 
surface  discharge  so  created  contrasts  with  the  volumetric 
effect  observed  when  the  electrodes  are  separated  by  the 
huid.  Typically,  the  actuator  is  excited  by  powering  the  ex¬ 
posed  electrode  at  radio  frequency  (rf)  voltages,  while  the 
embedded  electrode  is  grounded.  In  another  arrangement, 
both  electrodes  are  powered  with  signals  separated  by  a  beat 
frequency.  The  excitation  induces  a  complex  unsteady  inter¬ 
action  between  the  two  electrodes  and  the  huid,  details  of 
which  depend  on  frequency,  voltage,  geometric  conhgura¬ 
tion,  and  dielectric  constants  of  the  media.  In  the  surface 
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discharge,  within  a  very  short  time  after  breakdown,  the 
charge  buildup  at  the  dielectric  surface  sets  off  microdis¬ 
charges  of  nanosecond  duration,  limiting  the  electric  held  at 
the  location  of  the  microdischarge  such  that  the  charge  cur¬ 
rent  at  this  position  is  cut  off.  Experimental  evidence  shows 
that  there  is  no  runaway  state  for  the  parameters  under  con¬ 
sideration  and  that  an  asymptotic  (quasi)  periodic  state  is 
reached,  with  a  dominant  frequency  that  is  locked  to  the 
input  perturbation.  For  a  given  interelectrode  distance,  as  the 
applied  voltage  becomes  sufficiently  large,  the  dielectric  sur¬ 
face  adjacent  to  the  rf  electrode  produces  a  barrier  discharge, 
which  weakly  ionizes  the  surrounding  gas.  The  combination 
of  electrodynamic  body  force  and  collisional  processes, 
whose  detailed  mechanics  remains  a  matter  of  current  re¬ 
search,  ultimately  transfers  momentum  acquired  from  the 
electric  held  by  the  charged  particles  to  the  neutrals  which 
are  the  primary  species. 

The  advantages  of  dielectric  barrier-based  discharges  in¬ 
clude  an  absence  of  moving  parts,  rapid  on-off  features,  and 
the  ability  to  apply  body  forces  in  a  relatively  precise  manner 
by  deploying  advanced  electromagnetic  technology.  The  de¬ 
vices  are  thus  optimally  suited  for  how  control  in  wall  layers 
or  separated  layers1-10  The  main  parameters  applied  in  ex¬ 
perimental  observations  include  peak-to-peak  voltage  be¬ 
tween  2-20  kV  at  1-50  kHz  rf,  which  are  suitable  for  ac¬ 
tuation  at  atmospheric  pressure  at  low  speeds  to  0(10)  Torr 
at  high  speeds.  Specifically,  at  high  pressures  the  huid  is 
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FIG.  1.  Schematic  of  the  rf  induced  dielectric  barrier  discharge  device.  The  shaded  region  is  an  insulator  with  embedded  (grounded,  (f>=0)  and  exposed 
[powered,  (f>0  f (cot)]  electrodes  shown  in  thick  solid  line  segments,  (a)  Square  (S)  geometry  with  alumina  insulator,  (b)  Rectangular  (R)  geometry  with  Kapton 
insulator. 


highly  collisional,  inducing  an  efficient  momentum  transfer 
between  charged  and  neutral  species. 

The  collisional  rf  sheath  dynamics  of  near-surface  re¬ 
gions  is  fundamentally  different  from  that  under  direct  cur¬ 
rent  (dc)  and/or  collisionless  conditions  and  substantially 
more  difficult  to  simulate  because  of  their  unsteady  nature. 
While  dc-based  methods  are  a  useful  study  for  energy  inter- 
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actions  of  an  already  ionized  flow,  recent  literature 
shows  that  energy  budgets  will  depend  crucially  on  dynamic 
nonequilibrium  ionization  techniques.  The  load  factor  for  dc 
sheath  application  is  of  order  1,  far  from  the  most  energy 
efficient  Stoletow  point,  thus  generally  unsuitable  for  ioniza¬ 
tion  purposes.  While  the  earliest  atmospheric  glow  discharge 
research  by  von  Engle  et  al.u  was  published  over  70  years 
ago,  recent  experimental  successes  of  rf-induced  glow  dis¬ 
charge  using  air  at  atmospheric  condition1-6  have  triggered 
tremendous  interest  in  the  electro-gas  dynamic  flow  control 
and  energy  management  community.  Microfilaments  of 
nanosecond  duration  with  many  current  pulses  in  a  half  cycle 
maintain  the  optical  glow  in  dielectric  barrier  discharge 
(typically,  it  is  the  light  emission  spectrum  that  remains  dif¬ 
fused  to  create  a  uniform  glow).  The  small  time  scale  of  the 
discharge  limits  charge  transport  and  induces  minimal  heat¬ 
ing  of  the  bulk  gas.9,15  Thus,  most  of  the  electric  energy  is 
utilized  for  exciting  the  carrier  gas,  resulting  in  efficient  ap¬ 
plications  to  boundary  layer  flow  actuation.  Interestingly,  one 
atmosphere  uniform  glow  discharge  plasma  (OAUGDP)1,4 
also  shows  similar  flow-induction  phenomenology,  but  ap¬ 
parently  without  the  filamentation.  Henceforth,  we  will  use 
the  term  DBD  loosely  to  address  optical  glow  discharge  at 
high  pressures. 

The  typical  distance  between  dielectric  barrier  discharge 
(DBD)  electrodes  is  —0.1  mm  in  plasma  displays,  over 
1  mm  in  ozone  generators,  to  several  centimeters  in  C02 
lasers.15  It  has  been  observed  that  in  some  gases  (e.g.,  he¬ 
lium)  the  discharge  glow  is  stable,  whereas  for  the  same 
interelectrode  gap  and  pressure  air,  nitrogen,  oxygen,  and 
argon  transition  into  an  unstable  filamentary  discharge.15 
While  there  are  two  types  of  DBD,  namely  volume  discharge 
between  two  parallel  plates  and  surface  discharge  over  a  di¬ 


electric  boundary,  this  study  focuses  on  the  latter. 

Despite  the  extensive  recent  experimental  data  on  atmo¬ 
spheric  DBD  interactions  for  flow  actuation,  the  theoretical 
modeling  effort  is  still  in  its  infancy.  These  few  efforts  span 
a  range  of  phenomenological  to  rudimentary  first-principle- 
based  methods.4,5,7-10,16  For  example,  Roth4  utilized  a  postu¬ 
lated  balance  between  the  electrostatic  force  and  the  hydro- 
dynamic  force  to  estimate  the  induced  gas  velocity.  Enloe  et 
al.5  have  analytically  shown  in  a  simple  one-dimensional 
axisymmetric  configuration  that  two-dimensional  electric 
force  modeling  is  essential  to  the  understanding  of  the  glow 
discharge  flow  actuation  process.  For  accuracy  and  fidelity,  it 
is  imperative  that  the  force  model  be  derived  from  first  prin¬ 
ciples  through  a  simulation  of  the  elementary  mechanisms 
that  yield  the  discharge.  Shyy  et  al.7  developed  an  empirical 
model  based  on  pictures  of  the  glow.  Although  rudimentary, 
it  provides  a  spatial  variation  of  force  through  an  assumed 
triangular  distribution  of  the  electrostatic  body  force  down¬ 
stream  of  the  electrode.  Recently,  Gaitonde  et  al.16  param¬ 
etrized  the  force  direction,  and  incorporated  it  into  fully  3D 
direct  numerical  simulations  to  described  the  response  of  the 
flow  past  a  stalled  NACA  0015  airfoil  at  15°  angle  of  attack. 

For  accuracy  and  fidelity,  however,  it  is  imperative  that 
the  force  model  be  derived  from  first  principles  through  a 
simulation  of  the  elementary  mechanisms  that  yield  the  dis¬ 
charge  characteristics  and  flow  actuation  similar  to  that 
found  in  the  experiments.  Roy  and  Gaitonde  demonstrated 
such  a  model,  first  for  volume  discharge  between  two  dielec¬ 
tric  coated  electrodes8  and  then  for  capturing  effects  of  sur¬ 
face  discharge  on  quiescent  helium  gas,9,10  using  a  finite- 
element-based  multidimensional  multifluid  formulation  of 
plasma  sheath  at  atmospheric  conditions.  While  the  first 
model  was  relatively  basic,  the  two-dimensional  model  of 
high  pressure  surface  glow  discharge  documented  a  consis¬ 
tent  first-principles  formulation  of  charge  and  neutral  number 
densities,  their  momentum  dynamics,  electrostatic  field,  and 
potential  distribution  for  an  asymmetric  DBD  configuration. 
The  model  treated  the  insulator  and  the  gas  simultaneously 
and  integrated  the  Poisson  equation  into  the  charge  and  gas 
dynamics.  Simulations  showed  reasonable  comparison  with 
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reported  rf  glow  discharges17  in  partially  ionized  helium  gas 
between  two  insulated  electrodes.  The  current  effort,  focused 
as  it  is  on  surface  discharges  which  have  proven  more  attrac¬ 
tive  in  flow  control,  is  a  logical  next  step  in  the  evolution  of 
computational  methods  for  high  pressure  glow  discharges. 

The  present  paper  reports  the  first  two-dimensional 
tightly  coupled  first-principles  model  demonstrating  control 
of  a  separated  flow  at  atmospheric  pressure.  The  methodol¬ 
ogy  ensures  a  consistent  calculation  of  charge  and  neutral 
number  densities,  their  momentum  dynamics,  electric  field, 
and  potential  distribution.  The  formulation  is  anchored  in  an 
efficient,  robust,  and  versatile  module-based  finite-element 
code,  to  facilitate  the  prediction  of  the  electric  field  from  the 
charge  distribution  computed  from  the  grounded  and  applied 
rf  potential  and  the  dielectric  properties  of  insulator  and  gas 
media.  The  electric  force,  computed  as  a  product  of  self- 
consistent  charge  and  electric  field,  is  employed  to  calculate 
the  flow  actuation  effect.  For  simplicity,  here  again  helium  is 
used  as  a  carrier  gas  because  its  plasma  coefficients  and 
chemistry  are  known.  These  self-consistent  simulations  lay 
the  foundation  for  a  detailed  analysis  of  the  interaction  be¬ 
tween  charged  and  neutral  particles  in  the  important  situation 
where  separation  is  affected.  This  facilitates  a  more  insight¬ 
ful  understanding  of  the  complex  relation  between  the  body 
force,  inertial,  viscous,  and  pressure  fields  to  augment  ex¬ 
perimental  observations.  This  work  also  provides  a  basis 
through  which  future  efforts  can  investigate  the  impact  of 
dimensionless  parameter  variation,  particularly  of  the  Rey¬ 
nolds  number,  to  yield  guidance  on  the  potential  scalability 
of  these  devices. 


II.  MODEL  DETAILS 


The  following  two-dimensional  three- species  collisional 
plasma-sheath  model  is  solved: 


Charge  continuity: 


fag 

dt 


dnJJ, 


9H 


dX; 


=  nez  -  rneni} 


(la) 


Charge  momentum: 


d<p  dna 

fJXj  dXj 


(lb) 


(#<!>\ 

Potential:  e  — y  =e(ne-nl).  (2) 

\  axj ! 

The  charge  particle  a=e,i  distributions  are  considered  non- 
Maxwellian.  The  electron  temperature  is  nearly  uniform  at 
1  eV  =  1 1  600  K  and  the  ions  and  neutrals  are  in  thermal 
equilibrium  at  300  K.  The  working  gas,  helium,  is  main¬ 
tained  at  bulk  pressure  p  =  300  Torr  and  temperature  T 
=  300  K.  For  the  operating  condition  of  interest,  the  inertia 
terms  are  small  relative  to  the  effect  of  the  electrical  field  and 
collisional  interactions,  and  thus  neglected.  The  collisional 
momentum  exchange  is  vital  to  the  charge-to-neutral  mo¬ 
mentum  transfer. 

The  electron  diffusion  is  obtained  from  the  Einstein  re¬ 
lation,  De={Tel e)jxe,  where  Te  is  the  energy  in  electron  volts, 
e  is  the  elementary  charge,  s  is  the  permittivity,  and  jme 


=  e/(meveh)  is  mobility  of  an  electron,  where  veh~ven 
~  1012/s  is  the  electron-neutral  collision  frequency.  In  this 
context,  we  note  that  transport  properties,  which  have  been 
taken  from  the  literature,  are  to  be  viewed  as  nominal  values 
facilitating  the  development  of  the  numerical  framework, 
rather  than  to  fill  gaps  in  current  thermo-chemical  data  for 
the  environment  of  interest.  In  this  same  vein,  is  the  ion 
diffusion18  in  cm2/s  at  300  K,  and  the  ion  mobility  is 
given  as19 

p/uit  =  8  X  103(1  -  8  X  10 ~3E/p)  cm2  V-1  s-1  Torr, 
for  Elp  ^  25  V  cm-1  Torr-1, 


PVi  = 


4.1  X  104 

yEp 


27.44  \ 
(. Elp)312) 


cm2  V  1  s  1  Torr, 


for  E/p  >  25  V  cm  1  Torr  \ 


(3) 


where  E-  ^E^+Ey  is  the  electric  field  magnitude  with  com¬ 
ponents  Ej=-d(p/  dxj  and  p  is  the  gas  pressure.  The  index  j 
takes  the  values  v  and  y.  Note  that  the  ion  mobility  is  based 
on  an  empirical  relation  that  contains  the  electric  field  and 
pressure,  assuring  consideration  of  collisional  effects.  How¬ 
ever,  quantitative  results  will  depend  on  the  nominally  cho¬ 
sen  values.  The  ionization  rate  z  for  helium  gas  used  here 


z=Aexp((£^°r^jEs"1; 

A  =  4.4  cm-1  Torr-1  and  B-  14  [V/(cm  Torr)]0,4, 


(4) 


where  ne  is  the  electron  number  density.  The  coefficient  of 
two-body  and  three-body  recombination  on  the  right-hand 
side  of  Eq.  (la)  is  given  as20  r=[VethaT(Veth)]  cm3/s.  Here, 
Veth  is  the  electron  thermal  velocity  and  crei  is  the  electron- 
ion  collision  cross  section.  For  simplicity,  we  have  negligible 
secondary  emission  at  this  operating  frequency.  The  follow¬ 
ing  boundary  conditions  are  maintained:  At  the  electrode,  we 
maintain  De(dne/ dx2) =nepLe(d<fi/ dx2)  -  crne  Veto  and 

Di(dni/dx2)  =  0.  The  factor  cr  is  set  to  0  or  1  depending  on  the 
positive  or  negative  value  of  the  exposed  electrode  in  the  rf 
cycle. 

Importantly,  any  balance  between  the  electrostatic  and 
fluid  dynamic  pressure  gradients  cannot  be  guaranteed  for  a 
fundamentally  unsteady  process.  We  use  the  electrodynamic 
body  force  e(nrne)YL  (although  no  Lorentz  force  is  present, 
the  force  is  created  by  charge  transport)  to  compute  the  av¬ 
eraged  velocity  of  the  working  fluid  is  based  on  the  follow¬ 
ing  continuity  and  momentum  equations: 

dpldt+dpVfjldXj  =  0,  (5  a) 


?Xii 

dt 


+ v, 


.W£L 

1  dxj 


=  e(nt  -  ne)Ej  - 


RT  dp  d 

+  ri 

p  dxj  dxt 


(5b) 


where  Tjt= pi/ p(dV ^1  dxj+  dV fj/ dxj)-\{fi/ p)3V fk/ dxkSjh  with 
/x  as  the  gas  viscosity  (=2.066  X  10-5  N  s/m2)  and  Sjl  is  the 
Kronecker  delta.  Helium  is  assumed  as  an  ideal  gas.  The 
factor  uj  is  introduced  to  modulate  the  effect  of  electric  field. 
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The  multiscale  ionized  gas  (MIG)  flow  code  was  used  to 
solve  Eqs.  (la),  (lb),  (2)-(4),  (5a),  and  (5b).  MIG  has  a 
module-based  computational  platform  that  implements  the 
high-fidelity  finite-element  (FE)  procedure  adapted  from 
fluid  dynamics  to  overcome  the  stiffness  of  the  equations 
generated  by  multispecies  charge  separation  phenomena. 
This  code  has  been  used  and  validated  for  several  fluid 
dynamics,  ionized  gas,  ’  and  micro/nanoscale  ’  flow 
problems.  FE  techniques  are  especially  suitable  for  their 
adaptability  to  arbitrary  multidimensional  geometries  and 
boundary  conditions.  Here,  a  2D  biquadratic  FE  formulation 
is  employed  along  with  fourth-order  Runge-Kutta  time 
marching  to  solve  the  continuity  and  momentum  equations 
for  all  species  and  fluid.  The  solution  process  consists  of  two 
steps.  First,  using  Eqs.  (la),  (2),  and  (3),  a  global  matrix  is 
formed  and  solved  simultaneously,  obviating  the  need  for 
any  special  subiteration  for  the  Poisson  solver.  The  charge 
density  and  electric  field  thus  calculated  are  then  used  in  Eqs. 
(5a)  and  (5b)  to  predict  the  fluid  velocity.  The  Galerkin  weak 
statement  associated  with  a  variational  integral  underlines 
the  development  of  this  numerical  algorithm.  Note  that  the 
selection  of  a  test  function  orthogonal  to  the  trial  function  in 
Galerkin  formulations  guarantees  the  minimization  of 
error,25  making  it  better  suited  for  these  problems.  Details  of 
the  MIG  code  are  reported  elsewhere  (for  example,  see  Refs. 
20  and  21).  For  the  species  equations  (l)-(4),  a  solution 
adaptive  time  scale  of  up  to  100  times  the  dielectric  relax¬ 
ation  time  scale26  is  adopted  to  calculate  the  charge  separa¬ 
tion  and  electric  field  up  to  a  periodic  asymptote.  The  gas 
equations  (5  a)  and  (5b)  are  then  solved  using  the  volume 
specific  electrodynamic  body  force  to  predict  the  fluid  veloc¬ 
ity  and  continuity  evolving  in  time  with  the  flow  time  scale 
of  ~ms.  The  solution  of  the  Newton-Raphson  iteration  is 
converged  at  any  given  time  step  when  the  maximum  value 
of  the  residual,  relative  L2  norm  for  each  of  the  state  vari¬ 
ables,  becomes  smaller  than  a  chosen  convergence  criterion 
of  10“4. 

III.  PROBLEM  DESCRIPTION 

Two  distinct  geometries  denoted  as  S  (square)  and  R 
(rectangle)  are  simulated. 

A.  Geometry  S 

The  two-dimensional  computational  domain 
(x,y:  0,1  cm)  consists  of  a  lower  half  (y:  0,  0.5  cm)  insu¬ 
lator  with  dielectric  constant  ed,  and  a  upper  half  filled  with 
inert  helium  gas  of  Sf=  1.0055  e0,  where  e0  is  permittivity  of 
vacuum.  Inside  the  insulator  the  current  due  to  motion  of 
charged  particles  is  forced  to  zero,  while  the  displacement 
current  is  balanced  with  the  total  current  at  the  gas -dielectric 
boundary.  The  bottom  encapsulated  electrode  is  grounded 
and  a  rf  alternating  frequency  of  5  kHz  with  peak-to-peak 
potential  of  2  kV  is  imposed  at  the  top  electrode.  Each  elec¬ 
trode  is  infinitesimally  thin  and  2  mm  long.  In  order  to  iden¬ 
tify  the  effect  of  electrode  configuration  two  different  elec¬ 
trode  arrangements  were  simulated.  In  the  first,  denoted  as 
S-l,  the  two  electrodes  were  kept  vertically  aligned.  Hereaf¬ 
ter,  we  will  denote  this  as  the  symmetric  configuration.  In  the 


Stream  wise 


FIG.  2.  Representative  computational  mesh  shows  highly  nonuniform 
nature. 


second  case,  denoted  as  S-2,  the  electrodes  are  offset  hori¬ 
zontally  by  2  mm.  This  will  be  noted  as  the  asymmetric  con¬ 
figuration.  The  gas  is  assumed  initially  stagnant  (V^-= 0).  The 
solutions  are  verified  by  qualitative  comparison  with  the  re¬ 
ported  results.  The  results  are  then  employed  to  explore  the 
enhancement  of  near  wall  fluid  velocity. 

B.  Geometry  R 

The  two-dimensional  computational  domain 
(-1.252,3.880  cm)  X  (-0.1 ,2.0  cm)  consists  of  a  Kapton 
polyimide  insulator  in  the  lower  half  (y:  -0.1,0  cm)  with  a 
dielectric  constant  sd=  3.5  e0  and  the  upper  half  (y:  0,2  cm) 
filled  with  inert  helium  gas  of  e^=  1.0055  s0,  where  s0  is 
permittivity  of  vacuum.  Inside  the  insulator  the  current  due 
to  motion  of  charged  particles  is  forced  to  zero  while  the 
displacement  current  is  balanced  with  the  total  current  at  the 
gas-dielectric  boundary.  The  schematic  in  Fig.  1  shows  two 
electrodes  in  which  the  bottom  electrode  is  grounded  and  a  rf 
alternating  frequency  of  5  kHz  with  rms  potential  of  2  kV  is 
imposed  at  the  top  electrode.  Each  electrode  is  infinitesi¬ 
mally  thin  and  1.2  cm  long.  They  are  vertically  displaced  by 
1  mm  with  no  horizontal  overlap.  The  highly  stretched  rep¬ 
resentative  computational  mesh  shown  in  Fig.  2  consists  of 
51  X  53  two-dimensional  biquadratic  Lagrange  finite  ele¬ 
ments.  The  insulator  material  is  modeled  using  12X53  bi¬ 
quadratic  elements.  The  nodes  are  attracted  toward  the  wall 
and  electrode  edges  for  capturing  better  solution  details.  The 
highest  and  lowest  elemental  aspect  ratios  are  approximately 
90  and  10,  respectively.  In  order  to  identify  the  charge  char¬ 
acteristics  and  its  effect  on  gas  flow,  two  different  simula¬ 
tions  were  performed  on  a  single  geometric  configuration 
with  an  asymmetrically  displaced  electrode-pair  arrange¬ 
ment.  In  the  first,  no  gas  velocity  was  imposed,  i.e.,  quies¬ 
cent  flow.  Hereafter,  we  will  denote  this  as  case  R-l.  In  the 
second,  an  upstream  flow  of  lOm/s  was  imposed  at  v 
=  -1.252  cm  with  an  angle  of  attack  (AO A)  of  +12°.  This 
will  be  noted  as  case  R-2. 

For  all  simulations  (S-l,  S-2,  R-l,  and  R-2)  electrons  are 
assumed  to  be  isothermal  at  the  boundary  and  maintained  at 
1  eV  (—11  600  K)  while  the  ions  are  cold  (300  K)  at 
300  Torr.  The  electron  flux  at  the  exposed  electrode  is  based 
on  the  electron  thermal  velocity  and  is  neglected  unless  di¬ 
rected  toward  the  wall.  Homogeneous  Neumann  boundary 
conditions  are  applied  for  all  other  variables  at  other  bound¬ 
aries.  Zero  flux,  i.e.,  homogeneous  Neumann  boundary  con- 
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FIG.  3.  (Color  online)  (a),  (b)  Solution  for  symmetric  configuration,  (a)  Computed  axial  electric  field  contours  with  embedded  electric  field  lines,  (b) 
Streamwise  velocity  contours  calculated  from  the  electric  field  in  (a)  shows  counteracting  flow  field,  (c),  (d)  Solution  for  asymmetric  configuration,  (c) 
Computed  axial  electric  field  contours  with  embedded  electric  field  lines,  (d)  Streamwise  velocity  contours  calculated  from  the  electric  field  in  (c)  shows 
positive  flow  field. 


ditions,  are  applied  at  all  other  boundaries.  The  solutions  are 
verified  by  qualitative  comparison  with  the  reported  results. 
The  results  are  then  employed  to  explore  the  enhancement  of 
near-wall  fluid  velocity.  No-slip  gas  velocity  boundary  con¬ 
dition  is  imposed  for  all  cases.  Also,  the  insulator  is  consid¬ 
ered  impermeable  to  gas. 

IV.  RESULTS  AND  DISCUSSION 

Figure  3  displays  several  features  of  the  numerical  pre¬ 
dictions  on  both  the  symmetric  S-l  and  asymmetric  S-2  con¬ 
figurations.  For  S-l,  Fig.  3(a)  shows  the  axial  electric  field 
contours  with  embedded  electric  field  lines  computed  for  the 
symmetric  configuration  (at  ~7t/2  radians).  The  electric 
field  lines  trace  the  vector  path  moving  in  a  trajectory  from 
the  exposed  instantaneous  anode  to  the  grounded  cathode 
inside  the  dielectric  in  a  symmetric  fashion.  The  ions  travel 
along  these  repelling  field  lines  of  equal  strength  pushing 
outward  from  the  center  above  the  rf  electrode  inducing  the 
neutral  gas  particles  to  follow  suit.  The  net  effect  is  clearly 
seen  in  the  streamwise  gas  velocity  contours  in  Fig.  3(b) 
calculated  from  the  electric  field  depicting  a  counteracting 
flow  field  that  reduces  the  neutral  gas  density  in  the  region 
above  the  exposed  (rf)  electrode.  An  interesting  outcome  is 


that  the  fluid  over  the  electrode  will  respond  to  this  pulsating 
field  by  forming  a  jet,  whose  oscillations  will  depend  on  fluid 
inertia  and  the  instantaneous  polarity. 

However,  this  characteristic  behavior  changes  when  the 
embedded  electrode  is  displaced  axially  relative  to  the  ex¬ 
posed  electrode.  Figure  3(c)  plots  the  instantaneous  (at 
~7t/2  radians)  axial  electric  field  distribution  for  an  exposed 
anode  and  a  grounded  cathode  in  asymmetric  configuration. 
The  magnitude  of  the  axial  field  shows  a  large  round  peak  (in 
red)  at  the  right  of  the  anode  and  a  small  negative  (in  blue) 
just  left  of  the  anode.  The  electric  field  lines  overlaid  on  this 
clearly  represent  a  directional  bias  toward  the  right.  The 
dominant  acceleration  of  ions  in  this  direction  occurs  along 
these  lines,  leaving  only  a  relatively  smaller  regime  moving 
in  the  opposite  direction.  The  net  result  is  an  order  of  mag¬ 
nitude  higher  (than  the  symmetric  case)  gas  velocity  down¬ 
stream  of  the  right  edge  of  the  exposed  anode  with  a  two 
orders  of  magnitude  smaller  negative  velocity  near  the  left 
edge  [see  Fig.  3(d)].  These  solutions  agree  with  the  experi¬ 
mental  findings27  of  the  existence  of  a  low  pressure  region  in 
the  plasma  and  flow  induction  by  both  symmetric  and  asym¬ 
metric  plasma  actuators.  Also,  the  experiments  in  Ref.  27 
show  a  peak  velocity  1  mm  above  the  exposed  electrode 
which  matches  with  the  numerical  prediction.  Simulation  re- 
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FIG.  4.  (Color  online)  Time  evolution  of  computed  axial  electric  field  (volts/cm)  distribution  in  the  computational  domain. 


suits  highlight  that  the  asymmetry  enhances  the  electric  field 
and  thus  influences  flow  direction. 

A  clearer  understanding  of  the  asymmetric  configuration 
S-2  may  be  derived  by  exploring  various  quantities  at  several 
phases  of  the  excitation.  The  time  evolution  of  the  potential 
distribution,  axial  electric  field,  charge  distribution,  and 
stream  wise  and  crosswise  electron  current  components  are 
shown  in  Figs.  4-6  for  four  time  stations  inside  a  half  cycle. 
The  axial  electric  field  in  Fig.  4  shows  the  development  of  a 
near-circular  growth  of  positive  peak  of  approximately  1.1 
X  105  V/cm  and  negative  trough  of  5.5  X  104  V/cm  at  the 
right  and  left  end  of  the  exposed  electrode,  respectively.  The 
computed  solutions  are  very  similar  to  the  experimental  data 
of  Corke  et  al.  and  Roth  et  al. ,  where  the  exposed  elec¬ 
trode  is  situated  upstream  of  the  peak  location  of  the  electric 
field  and  thus  the  peak  of  velocity.  Note  that  the  highest 
value  of  the  electric  field  about  the  edges  of  the  grounded 
electrode  is  at  least  twice  as  much  of  that  near  the  exposed 
electrode.  The  peak  of  the  crosswise  component  Ey  (not 
shown  here  due  to  space  limitations)  is,  however,  just  above 
the  electrodes.  Specifically,  above  the  positive  electrode  Ey 
takes  values  of  about  4  X  104  V/cm. 

Based  on  the  Poisson  equation  (2),  the  divergence  of  the 
electric  field  is  directly  related  to  the  charge  distribution.  The 
nonuniform  nature  of  the  distribution  along  the  edge  of  the 
dielectric  is  further  elucidated  in  the  two-dimensional  charge 
solutions  plotted  in  Fig.  5.  Note  that  a  residual  negative 
charge  remains  deposited  on  the  dielectric  surface  down¬ 


stream  of  the  exposed  electrode  for  the  entire  duration  of  the 
cycle.  This  causes  a  net  electric  force  in  the  direction  from 
the  electrode  to  the  downstream  surface.  This  is  the  funda¬ 
mental  reason  why  the  net  near- wall  jet  velocity  is  unidirec¬ 
tional.  The  negative  charge  accumulation  at  the  surface  is 
consistent  with  the  physical  model  of  the  DBD  for  the  charge 
deposition  at  the  surface  of  the  dielectric.5 

Figure  6  displays  the  peak  streamwise  electron  current 
(+2.4, -4.1  mA)  observed  at  the  two  corners  (right  and  left) 
of  the  exposed  electrode  with  a  peak  value  of  -28.1  mA  at  i t 
radians,  and  the  minimum  current  (+0.3, -0.1  mA)  is  com¬ 
puted  at  it  1 2  radians  of  the  phase.  An  explanation  for  this 
behavior  may  be  proposed  based  on  arguments  pertaining  to 
the  observed  self-limiting  nature  of  the  discharge.  In  surface 
DBD  there  is  no  real  current  path  between  the  electrodes. 
However,  the  charge  buildup  on  the  insulator  surface  op¬ 
poses  the  applied  voltage,  shutting  off  the  discharge  just  be¬ 
yond  the  peak  when  the  applied  voltage  starts  decreasing, 
thus  self-limiting  the  system.  This  is  also  confirmed  by  the 
experimental  data  (see,  for  example,  Fig.  2  of  Ref.  17), 
where  the  current  magnitude  peaks  at  about  tt  radians.  The 
local  fluctuations  above  the  rf  electrode  as  seen  in  Fig.  6(b) 
show  a  possible  electrical  double  layer  about  the  sheath 
edge.  The  crosswise  electron  current  in  Fig.  8  also  shows 
similar  surface  characteristic  with  low  currents  at  tv  1 2  radian 
and  large  peak  value  (-243  mA)  near  tt  radian. 

Figure  7(a)  plots  the  streamwise  component  of  the  com- 
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FIG.  5.  (Color  online)  Time  evolution  of  charge  (nf-ne)  distribution  in  the  computational  domain. 


puted  gas  velocity  depicting  a  strong  peak  downstream  of  the 
rf  electrode  away  from  the  dielectric  surface.  The  initial  con¬ 
dition  is  zero  flow,  which  makes  the  computational  problem 
more  challenging.  The  local  vertical  line  plots  at  seven 
streamwise  locations  describe  how  the  flow  velocity  in¬ 
creases  in  response  to  the  increase  in  axial  electric  field  (Fig. 
3).  The  peak  of  the  wall  jet  like  feature  occurs  at  —2  mm 
downstream  of  the  exposed  electrode,  beyond  which  diffu¬ 
sive  momentum  transfer  increases  the  height  of  the  region 
influenced  by  the  DBD,  and  the  peak  velocity  diminishes. 
The  velocity  induced  by  the  discharge  is  relatively  small,  but 
should  increase  if  the  interelectrode  gap  is  decreased  in  both 
vertical  and  streamwise  directions  as  that  will  increase  the 
magnitude  of  the  force  generated  by  the  electric  field.  In 
Refs.  2  and  4,  it  has  been  conjectured  that  it  is  the  Lorentzian 
collision  between  the  ions  and  neutrals  that  brings  about  the 
plasma-to-gas  momentum  transfer.  Our  results,  in  fact,  indi¬ 
cate  that  for  helium  the  electron-neutral  collisional  influence 
meneven(Vej~Vnj)  1 n  this  process  may  not  be  negligible  as 
compared  to  the  effect  of  ion-neutral  collision,  MniV^iYij 
-Vnj).  It  is  because  M>me  and  niVinV ij<nevenV ej  results  in 

mene  ^en  ^ e j  ~  Mil i  Vjn  Vf j . 

For  future  design  purposes,  it  is  important  to  understand 
the  effect  of  permittivity  of  the  dielectric  material  on  the 
electric  field  and  thus  on  the  flow  control.  Figure  8  shows  the 
effect  for  two  different  dielectric  constants  that  are  orders  of 
magnitude  apart  (although  realistic  materials  have  dielectric 
constants  below  10,  we  intend  to  amplify  the  effect).  The 
encircled  region  about  the  insulator  surface  highlights  the 


observation  that,  as  the  dielectric  constant  increases,  the 
crosswise  electric  field  decreases  noticeably,  bending  the 
electric  field  lines  almost  parallel  to  the  surface.  This  in¬ 
creases  the  current  component,  due  to  higher  resistivity  in 
the  electric  path  resulting  in  higher  power  requirement.  More 
investigations  are  required  to  confirm  this  inference. 

The  analyses  of  R-l  and  R-2  are  reported  in  Figs.  9-12. 
In  general,  the  electrohydrodynamic  (EHD)  body  force  cal¬ 
culated  based  on  the  product  of  charge  and  electric  field 
shows  very  strong  forward  and  downward  moving  vectors 
downstream  of  the  exposed  electrode  for  the  positive  part  of 
the  cycle.  At  the  negative  peak  of  the  cycle,  the  force  vectors 
concentrate  just  upstream  of  the  right  edge  of  the  powered 
electrode  with  majority  forces  going  backward  and  down¬ 
ward.  However,  for  the  rest  of  the  cycle  the  force  is  much 
smaller  as  compared  to  the  positive  part  (—10%,  42%,  and 
2%  of  the  peak  magnitude,  respectively,  at  7 r,  37t/2,  and  277 
radians).  The  average  body  force  will  thus  push  the  bulk  gas 
in  the  forward  and  downward  directions.  Recent  numerical 
simulation6,7  documented  the  gas  velocity  downstream  of  the 
right  edge  of  the  exposed  anode  with  a  two  orders  of  mag¬ 
nitude  smaller  negative  velocity  near  the  left  edge  of  the 
exposed  electrode  similar  to  what  we  observe  here. 

Figure  9  plots  the  streamwise  component  of  the  time 
average  of  volume  specific  body  force  for  case  R-l.  The  line 
trace  of  the  force  vectors  is  showing  a  directional  bifurcation 
(similar  to  that  reported  in  a  recent  PIV  experiment,28  Fig.  2) 
just  downstream  of  the  exposed  electrode.  This  will  tend  to 
guide  the  flow  in  both  directions  at  just  about  x=2cm  loca- 
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FIG.  6.  (Color  online)  Time  evolution  of  streamwise  component  of  the  electron  current  in  mA. 


tion.  The  positive  peak  very  close  to  the  dielectric  surface  is 
three  times  higher  than  the  negative  peak.  Due  to  fluid  iner¬ 
tia,  the  bulk  gas  will  only  respond  to  this  average  force, 
which  will  ensure  its  net  forward  motion.  The  momentum 
thus  imparted  to  the  gas  will  induce  a  velocity  along  the 
dielectric  surface.  Figure  10  plots  the  streamwise  component 
of  the  gas  velocity  computed  from  Eq.  (4)  for  case  R-l  at  six 
local  vertical  line  plots  upstream  and  downstream  of  the 
electrode  edge  and  shows  a  wall  jet  like  feature.  The  zero 
flow  initial  condition  makes  the  computational  problem  more 
challenging.  The  peak  of  this  wall  jet  of  2.4  m/s  occurs  at 
6  mm  downstream  of  the  exposed  electrode,  beyond  which 
diffusive  momentum  transfer  increases  the  height  of  the  re¬ 
gion  influenced  by  the  DBD,  and  the  peak  velocity  dimin¬ 
ishes.  Similar  velocity  profiles  have  been  observed  by  sev- 
eral  experimental  groups.  ’  ’  Note  that  for  helium  gas  at 
300  Torr  the  mobility  is  much  higher  than  atmospheric  air 
and  hence  the  predicted  induced  velocity  is  higher. 

Specific  details  of  gas  velocity  components  and  the 
streamwise  force  is  presented  in  Fig.  11  at  five  locations 
along  the  dielectric  at  t-  6  ms.  It  is  interesting  to  note  that 
the  peak  force  is  located  at  v=0.7  cm  in  the  middle  of  the 
initial  bubble,  while  the  peak  of  the  wall  jet  is  occurring 
downstream  of  the  exposed  electrode  at  x=2.0  cm.  The  lo¬ 
cation  of  the  actuator  is  thus  important  to  control  the  flow 
effectively. 

For  the  actuator  model  to  work  for  design  purposes,  it  is 
essential  to  test  it  for  low  speed  flow  conditions  where  a 
reasonable  amount  of  supporting  experiments  have  been  re¬ 


ported.  Figure  12  shows  the  transient  effect  of  DBD  actuator 
for  realistic  gas  flow  control.  The  time  scale  of  gas  flow  for 
the  present  dimensions  is  in  milliseconds  while  the  plasma 
time  scale  is  several  orders  of  magnitude  less.  For  this  case 
R-2  simulation,  first  an  initial  flow  condition  is  generated 
over  the  flat  plate  kept  at  a  12  deg  angle  with  respect  to  the 
helium  gas  inflow  without  the  rf  power  turned  on.  Figure 
12(a)  plots  the  stream  traces  of  gas  particles  based  on  the 
velocity  vectors  of  the  initial  flow  field  and  shows  a  strong 
separated  flow  (see  the  inset)  just  downstream  of  the  leading 
edge  of  the  plate.  Then,  the  power  is  turned  on  for  the  same 
electrode  arrangement  as  presented  earlier.  Figures 
12(b)— 12(d)  show  the  gradual  removal  of  the  recirculation 
bubble  due  to  the  separation  as  a  function  of  time.  The  strong 
positive  force  contours  (in  red)  shows  the  time  average  dy¬ 
namics  of  the  streamwise  component  as  it  literally  pushes  the 
recirculation  off  the  plate  within  6  ms.  The  fluid  reacts  to  the 
generated  mean  body  force  but  there  may  be  a  slow  interac¬ 
tion  with  the  local  plasma  generation  as  time  progresses.  To 
our  knowledge,  such  description  has  not  been  shown  in  the 
reported  literature  to  date. 

V.  CONCLUSIONS 

A  two-dimensional  finite-element-based  formulation  of 
plasma-fluid  interactions  is  given  for  a  partially  ionized 
plasma  using  the  multicomponent  fluid  equation.  The  model 
is  applied  to  simulate  an  atmospheric  surface  dielectric  bar¬ 
rier  discharge  for  partially  ionized  helium  gas.  The  computed 
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FIG.  7.  (Color  online)  (a)  Contours  of  the  streamwise  velocity  identify  the 
location  of  highest  momentum  transfer,  (b)  Calculated  streamwise  compo¬ 
nent  of  the  gas  velocity  distribution  at  different  locations  along  the  flow 
showing  similar  wall  jet  as  found  in  experiments  (Refs.  3,  4,  and  28). 
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results  are  similar  to  the  experimental  data  showing  that  the 
exposed  electrode  is  situated  upstream  of  the  peak  location 
of  the  electric  field  and  thus  the  imparted  momentum.  The 
two-dimensional  predictions  mimic  the  self-limiting  dis¬ 
charge  mechanism  due  to  the  charge  accumulation  on  the 
insulator  surface  and  the  lowest  electron  current  just  beyond 
the  peak  applied  voltage.  The  results,  however,  stress  that  the 
asymmetry  enhances  the  electric  field.  The  effect  of  increas¬ 
ing  insulator  dielectric  constant  in  the  computed  electric  field 
lines  shows  noticeable  decrease  in  the  crosswise  electric 
field,  possibly  increasing  the  power  requirement.  More  in¬ 
vestigation  is  needed  to  ascertain  optimum  electrical  and 
geometric  design  parameters.  In  our  simulations,  the  electric 
field  varies,  but  reaches  a  quasiperiodic  asymptote.  The  fun¬ 
damental  reason  for  the  net  near- wall  unidirectional  jet  ve¬ 
locity  can  be  explained  as  electrostatic  force  along  the  elec¬ 
tric  field  toward  the  electron  deposition  downstream.  The 
issue  of  whether  the  fluid  reacts  to  this  force  in  a  mean  or 
time-varying  sense  has  not  yet  been  resolved,  and  requires 
future  work.  Also,  we  have  utilized  300  Torr  pressure  while 
the  atmospheric  pressure  is  nearly  2.6  times  higher.  This  will 
increase  the  collisional  effects,  reducing  charge  mobility  and 
decreasing  resulting  induced  momentum  transfer  to  the  gas. 


FIG.  8.  (Color  online)  The  electric  field  lines  show  as  the  dielectric  constant 
increases  the  field  lines  at  the  surface  of  the  dielectric  bend  in  the  high¬ 
lighted  ellipse,  showing  possible  increase  in  current,  hence  the  power 
requirement. 

In  the  near  future,  the  model  will  be  extended  to  atmospheric 
air  with  negative  ions  and  additional  new  mechanisms  in  the 
source  terms.  A  model  for  realistic,  fully  three-dimensional 
geometric  and  electrode  configuration  is  also  in  develop- 


Streamwise  (cm) 

FIG.  9.  Time  average  of  streamwise  component  of  the  julN  force  per  unit 
volume  about  the  surface  of  the  actuator  shows  the  dominance  of  the 
streamwise  forward  (positive)  force  component.  Line  traces  show  the  direc¬ 
tion  of  force  which  will  guide  the  local  flow. 
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FIG.  10.  Computed  streamwise  velocity  induced  in  a  quiescent  helium  gas. 


ment,  as  is  an  exploration  of  the  effect  of  different  voltage 
shapes  (e.g.,  sawtooth,  square  wave,  etc.).  The  present  effort 
thus  provides  a  practical  tool  to  augment  experimental  obser¬ 
vations  in  exploring  flow  control  concepts  and  in  developing 
suitable  inputs  for  traditional  fluid  dynamics  codes  based  on 
the  Navier- Stokes  equations. 

A.  Major  future  advancements  needed 

The  electron  temperature  was  assumed  to  be  uniform  at 
Te=  1  eV.  This  assumption  can  lead  to  substantial  errors.  In 
fact,  at  these  high  pressures  (—300  Torr)  the  characteristic 
relaxation  time  of  electron  temperature  is  on  the  order  of  3 
nanoseconds;  accordingly,  the  relaxation  length  is  also  very 
short.  Therefore,  Te  is  a  function  of  the  local  value  of  E/N. 
Since  E  is  quite  nonuniform,  Te  is  also  quite  nonuniform. 
One  immediate  consequence  is  that  considerable  gradients  of 


FIG.  11.  Velocity  components  and  streamwise  volume  specific  force  distri¬ 
bution  at  specific  locations  (x  in  cm)  normal  to  the  dielectric  surface  after 
6  ms. 


Te  should  result  in  additional  electron  fluxes.  In  Eq.  (lb),  the 
diffusion  term  only  has  the  gradient  of  electron  density. 
However,  the  diffusion  term  should  in  general  be  propor¬ 
tional  to  the  gradient  of  electron  pressure  and  can  be  broken 
into  two  terms — one  with  electron  density  gradient,  and  the 
other  with  the  electron  temperature  gradient. 

(2)  Also,  electron  adsorption  on  the  dielectric  surface 
occurs  for  a  dielectric  coated  anode.  The  electron  avalanches 
flood  the  dielectric  surface  until  the  threshold  intensity  is 
attained  for  formation  of  streamers.  For  applied  rf  voltage, 
thermo-desorption  and  interaction  with  vibrationally  excited 
and  metastable  gas  molecules  may  cause  the  electron  emis¬ 
sion  from  the  dielectric  and  electrode  surfaces.  Maslov 


FIG.  12.  (Color  online)  Separation  control  for  an  incoming  flow  with  +12°  AOA  over  a  flat  plate  shows  the  gradual  flow  attachment  eradicating  the  separation 
bubble  in  a  few  milliseconds.  The  velocity  stream  traces  are  plotted  over  the  streamwise  specific  force  contours  showing  a  strong  average  force  over  the 
exposed  electrode. 
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et  al.  conjectured  that  the  secondary  electron  emission  due 
to  ion  current  seems  to  be  less  important  than  the  aforemen¬ 
tioned  mechanisms.  However,  the  electrons  emitted  from  the 
surface  become  the  source  for  electron  avalanche  develop¬ 
ment  and  their  further  transformation  into  streamers.  Thus 
for  completeness,  the  modality  of  plasma- wall  interaction 
should  be  carefully  investigated. 
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